Discrete Solitons and Vortices on Anisotropic Lattices 
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We consider effects of anisotropy on solitons of various types in two-dimensional nonlinear lattices, 
using the discrete nonlinear Schrodinger equation as a paradigm model. For fundamental solitons, 
we develop a variational approximation, which predicts that broad quasi-continuum solitons are 
unstable, while their strongly anisotropic counterparts are stable. By means of numerical methods, 
it is found that, in the general case, the fundamental solitons and simplest on-site-centered vortex 
solitons ("vortex crosses") feature enhanced or reduced stability areas, depending on the strength 
of the anisotropy. More surprising is the effect of anisotropy on the so-called "super-symmetric" 
intersite-centered vortices ("vortex squares"), with the topological charge S equal to the square's 
size M : we predict in an analytical form by means of the Lyapunov-Schmidt theory, and confirm 
by numerical results, that arbitrarily weak anisotropy results in dramatic changes in the stability 
and dynamics in comparison with the degenerate, in this case, isotropic limit. 



I. INTRODUCTION AND THE MODEL 

In the last two decades, nonlinear lattice (spatially dis- 
crete) systems have been a very rapidly growing area of 
interest for a variety of applications Such systems 
arise in physical contexts encompassing, inter alia, beam 
dynamics in coupled waveguide arrays in nonlinear op- 
tics j^, the time evolution of fragmented Bose- Einstein 
condensates (BECs) trapped in optical lattices (OLs) 0, 
coupled cantilever systems in nano- mechanics |j| , denat- 
uration of the DNA double strand in biophysics and 
even stellar dynamics in astrophysics [H . 

One of the main objectives of the research in this field 
is to achieve an understanding of intrinsically localized 
states (discrete solitons). In two-dimensional (2D) lat- 
tices, these are fundamental discrete solitons ,7] and dis- 
crete vortices (i.e., localized states with an embedded 
nonzero phase circulation over a closed lattice contour) 
H, ■ Most recently, a substantial effort was dedicated 
to the experimental creation of both these entities in 
photonic lattices induced in photorefractive crystals (al- 
though these systems are only quasi-discrete). In par- 
ticular, fundamental and dipole solitons, soliton trains 
and necklaces, and vector solitons have been reported 
[TQi] , as well as vortex solitons • Parallel developments 
in the experimental studies of soliton patterns in BECs 
have also been very substantial, leading to the creation of 
quasi-lD dark 12], bright 13] and gap 0] solitons. The 
generation of 2D BEG solitons in OLs has been theoreti- 
cally demonstrated to be feasible with the currently 
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available experimental technology j^. 

A paradigm dynamical lattice model that appears in 
the above-mentioned physical problems is the discrete 
nonlinear Schrodinger (DNLS) equation. Various ap- 
plications of the DNLS equation are well documented 

13, 3 • Besides being a generic asymptotic form of a 
whole class of lattice models (for small-amplitude nonlin- 
ear excitations) , it finds direct applications (where it fur- 
nishes extremely accurate description of the underlying 
physics) in terms of arrayed (ID) or bunched (2D) non- 
linear optical waveguides, BEGs trapped in strong OLs, 
and crystals built of optical or exciton microcavities. 

An interesting issue in this framework that has not 
received sufficient attention concerns the influence of 
anisotropy on the soliton dynamics in 2D lattices. 
Some of the settings mentioned above are inherently 
anisotropic, e.g., photorefractive crystals 0, 0, while 
others (in particular, the fragmented BEGs trapped in 
strong OLs 3]) can be easily rendered anisotropic by 
slight variations of control parameters, such as intensities 
of laser beams that create two sublattices which together 
constitute the 2D optical lattice. 

The aim of this paper is to understand how the lat- 
tice anisotropy affects 2D discrete solitons in the DNLS 
equation. Some findings reported below are surprising, 
demonstrating that anisotropy effects are not straight- 
forward. The straightforward expectation might be that 
weak anisotropy is a small perturbation that possibly al- 
ters details of parametric dependences of the observed 
phenomenology but does not change it "structurally" 
(i.e., essentially the same dynamical features as in the 
isotropic lattice occur, but at different positions in the 
parameter space). We find that for the simplest soli- 
ton and vortex structures this is indeed the case, while 
for more sophisticated ones it is not. More specifically. 
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we find that for especially symmetric (so-called "super- 
symmetric") vortices, with their center set at an inter- 
site position, and the topological charge equal to the size 
of the vortex square frame (see below for details), the 
isotropic lattice is a degenerate one, therefore even very 
weak anisotropy fundamentally alters the stability and 
dynamical properties of such structures. On the other 
hand, despite the delicate organization of the supersym- 
metric vortices, they constitute a structurally stable, i.e., 
physically meaningful, class of objects. 

We take the DNLS equation in the following form: 



(1) 



where Un,m{t) is the complex, 2D lattice field (the overdot 
stands for its time derivative), e is the lattice coupling 
constant, and 



Q! {Un+l,m + Un-l,m) + Un,m+1 
+Un,m-1 - 2(1 + a)u 



(2) 



is the anisotropic discrete Laplacian, which becomes 
isotropic with a = 1. Note that, unlike the continuum 
limit, no scaling transformation can cast the anisotropic 
DNLS equation into the isotropic form. Equation 
conserves two dynamical invariants: the Hamiltonian, 
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and norm, 



(4) 



where A is the frequency of the internal mode (equiva- 
lently the chemical potential in the context of BECs or 
the propagation constant in the context of optical waveg- 
uide arrays). 

Stationary solutions to Eq. ^ will be sought as 

Un.m = Mi°Lexp(iAi), (5) 
which leads to a stationary finite-difference equation. 
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(generally speaking, the discrete functions Mn°L may be 
complex). In the case of fundamental-soliton solutions, 
we will apply the variational approximation (VA) to the 
real version of Eq. (O, which is based on the fact that it 
can be derived from the Lagrangian, 
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After analyzing fundamental solitons by means of the 
VA, we will construct discrete solitons in the anisotropic 
model and will study their stability by means of numer- 
ical methods. For the numerical procedure, our start- 
ing point is always the anti-continuum (AC) limit corre- 
sponding to e = where configurations of interest 
can be constructed at will as appropriate combinations 
of on-site states, which are either Un,m = •\/Aexp(iAt) 
with A > at excited sites, and Un,m = at non-excited 
ones, cf. Eqs. Q and JHJl for the general case, e > 0. The 
stability of the solitons is then analyzed by linearizing 
Eq. (^) for perturbations around a stationary solution 



,(0) 



iAt 



(8) 



where S is an infinitesimal perturbation amplitude of the 
perturbation, and A is its eigenvalue. The Hamiltonian 
nature of the system dictates that if A is an eigenvalue, 
then so also are —A, A* and —A* (in the stable case, A is 
imaginary, hence this symmetry yields only two different 
eigenvalues, A and —A). Clearly, the stationary solution 
is unstable if at least one pair of eigenvalues features 
nonvanishing real parts. 

It is noteworthy that the instability against perturba- 
tions corresponding to purely real eigenvalues A in Eq. © 
can be predicted by the Vakhitov-Kolokolov (VK) crite- 
rion |l8j| : a soliton family, characterized by the depen- 
dence N{A) [recall N is the solution's norm defined by 
Eq. I0J], may be stable under the condition dN/dA > 0, 
and is definitely unstable in the opposite case. In par- 
ticular, this criterion (as well as the VA) was found to 
be very useful and quite reliable in the investigation of 
2D solitons in the Gross-Pitaevskii equation for BECs in 
2D and quasi- ID periodic OL potentials ^.nd even in 
2D quasi-periodic potentials (such as the Penrose tiling 
among others) 20]. 

Our study of different states in the anisotropic model 
and their properties is structured as follows. In Sec- 
tion II, we present the VA for fundamental solitons. In 
Section III, discrete solitons and vortex crosses with the 
topological charge S = 1 are considered, which are only 
perturbatively (weakly) affected by the anisotropy. In 
the following two sections, we will define and consider 
special "super-symmetric" configurations, with S = 1 
and S — 2, respectively, and compare them with simpler 
cases. Finally, in section VI we summarize findings and 
present our conclusions. 



II. VARIATIONAL APPROXIMATION FOR 
FUNDAMENTAL SOLITONS 

As was shown in Ref. for the one-dimensional 
DNLS equation (see also Ref. |2^, for a more rigorous 
variational approach applied to higher dimensional soli- 
tons in the isotropic case), the only analytically tractable 
variational ansatz for stationary fundamental solitons 
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may be based on the following cusp-shaped expression 
(in the 2D case, it has the shape of a cross cusp), 



,(0) 



Aexp (— a|n| — b\m\) . 



(9) 



with positive parameters a and b, that determine the 
widths of the soliton in the horizontal and vertical di- 
rections, and an arbitrary amplitude A. Note that ex- 
pression is indeed an exact solution to the linearized 
version of Eq. ©, which describes soliton tails, if A is 
linked to a and b by the dispersion relation. 



A = 2e [asinh2(a/2) +sinh2(5/2)] 



(10) 



The substitution of ansatz ^ makes it possible to cal- 
culate the corresponding effective Lagrangian explicitly. 
First of all, it is convenient to eliminate the amplitude 
in favor of the norm |0J. Indeed, the substitution of the 
ansatz in the definition of N yields ^ N tanh a-tanh b. 
After this, the effective Lagrangian becomes 
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Variational equations for the stationary profile are ob- 
tained from here in the form 



9Leff dLcs dL. 
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da 



db 



0. 



(12) 



In the general case, the explicit form of these equations 
is quite cumbersome (this will be treated numerically, 
see below). A detailed analysis is possible in two special 
cases, as specified below. 

First is the case of small a and (a, 6 <C 1), which im- 
plies broad solitons. Then, the expansion of the effective 
Lagrangian 111() yields 



L. 



off 



A 1 

-2^^+2^ 



ab+ -a'^b+ -ab"^ 

16e I 3 3 



5 ,A 5a A 
12^+12^ 



(13) 



and the variational equations H12|l following from Eq. H13|) 
generate the following solution: 
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(14) 
(15) 



As follows from these expressions, the underlying as- 
sumptions a, 6 ^ 1 indeed hold (i.e., the approximation 
is self-consistent) under the condition 



A« 



ae, if a ^ 1, 
e, if a > 1. 



(16) 



The broad (quasi-continuum) solitons predicted in this 
approximation are unstable according to the VK crite- 
rion, as Eq. 114|) immediately shows that dN/dA < 0. 

Note that the expansion of the dispersion relation H1U|) 
for the same case of small a and b yields aa^ + b^ — 
A/e. It is noteworthy that this relation, although derived 
independently of the variational equations, is consistent 
with Eq. ((T^ . 

Another tractable case is that of a strongly anisotropic 
soliton, which is broad (quasi-continuum) in either di- 
rection and narrow in the other, i.e., it corresponds to 
a ^ 1,6 ^ 1, or vice versa. If a is small and b is large, 
the variational equations H12|l yield the following results: 



, sinh I — 

3ae' V2 



-eaA. 



(17) 



These results are consistent with the underlying assump- 
tions (a ^ 1, 6 S> 1) under the conditions 



1< A/e < a. 



(18) 



On the contrary to the broad solitons given above by Eqs. 
(|14|l and H15|) . Eqs. H17|) show that the anisotropic solitons 
arc stable as per the VK criterion, as they obviously meet 
the condition dN/dA > 0. 

For the opposite strongly anisotropic case, with a ^ 1 
and 6 <C 1, the result is 



= -eA, 



(19) 



cf. Eqs. (|17|l . These expressions comply with the under- 
lying assumptions a 3> 1 , 6 <C 1 provided that 



a < A/e < 1, 



(20) 



cf . Eq. (fT^ . Similarly to the solution of Eq. ifTTjl , the one 
of Eq. H19|) obviously meets the VK stability criterion. 

Lastly, inequalities H18|) and (|20|) imply that the 
above solutions indeed pertain to the strongly anisotropic 
model, as the corresponding parameter a is large in the 
former case, and small in the latter one. We also no- 
tice that the condition H18() in the case of large a, or its 
counterpart H20|l in the opposite case of small a, is in- 
compatible with the respective condition (|16|l . i.e. (as 
one would expect), the existence regions of the unstable 
quasi-continuum solitons and stable strongly anisotropic 
ones have no overlap. 

For general a and b, the variational equations H12|l . 
with the effective Lagrangian (|ll|l . cannot be solved ex- 
plicitly and one has to find {N, a, b) solutions numerically 
for each (e. A) pair. In Fig. ^ we compare the results 
obtained from the VA with solutions obtained through 
numerically solving the stationary equation |(SJl. Fig. 
^a depicts the norm of the soliton solutions as a func- 
tion of the propagation constant A for several values of 
the anisotropy parameter a and for constant coupling 
(e = 1). As may be noticed from the figure, the VA (thin 



4 




0.5 .1 1.5 




0.5 1 £ 1.5 2 



FIG. 1: a) Norm of the solution vs. A for several values 
of the anisotropy and fixed coupling strength e — 1. For 
all the panels in this figure the anisotropy values are a = 
1.5, 1.25, 1, 0.75, respectively, for each curve from top to bot- 
tom. Thick lines (solid and dashed) represent direct numerical 
results and the thin lines represent the VA. The dashed lines 
correspond to unstable soliton solutions Note that the sign of 
the slope of A(A) reflects the stability of the soliton solutions 
as predicted by the VK criterion, b) The norm of the soliton 
solution as a function of the coupling strength for fixed A = 1. 
Once again, thick lines represent direct numerical results and 
thin lines illustrate the VA. 

lines) provides a good approximation to the actual solu- 
tion (thick lines). We also checked the stability of the 
constructed solutions by following the largest real eigen- 
value (see also details below), of the linearized problem 
defined in Eq. |(SJ). Stable solutions are depicted with 
solid lines while unstable solutions correspond to dashed 
lines. As is clear from the figure, the slope of iV(A) pre- 
dicts the stability of the solution according to the VK 
criterion (see above). Furthermore, since the VA gives 
a good approximation of N{K), it is possible to obtain 
a good estimate for the transition from stable to unsta- 
ble solutions, as A is decreased, using the VA together 
with the VK criterion. Finally, in Fig. QJb we fix A = 1 
and perform a similar calculation by varying the coupling 
strength e. Again, the VA (thin lines) approximates re- 
markably well the norm of the solutions (thick lines). 



III. FUNDAMENTAL SOLITONS AND VORTEX 
CROSSES: NUMERICAL RESULTS 

We start numerical computations with a single excited 
site in the AC limit, and continue the solution in e (for a 
fixed value of the anisotropy parameter a). The objective 



is to construct regular site-centered discrete solitons, with 
the anticipation that, as is known for the isotropic model 
(a = 1), the solitons will be stable up to a critical value 
of the coupling constant, i.e., at e < ecr [H, ll^. At 
e > Ccr, the discrete solitons are found to be unstable due 
to a real eigenvalue arising in the linearization around 
the soliton. In the numerical part of the work (unlike 
the VA considered above), we fix A = 1 in Eq. ||HJ), using 
the scaling invariance of Eq. and examine how e„ 
is affected by the variation of a. The results will be 
summarized in the form of two-parameter diagrams that 
chart regions of stable and unstable discrete states. 

For regular discrete solitons, such a diagram is pre- 
sented in Fig. 121 The top panel illustrates the fact 
that the increase of a gradually destabilizes the soli- 
tons, i.e., Ccr decreases with increasing a. Interestingly, 
the respective dependence is very well approximated by 
an empirical relation ecr = ^l\foL. More accurately, 
the best fit to this numerical dependence is given by: 
Ecr ~ 0.999q;~° ''8^. The middle panel in Fig.Elillustrates 
in more detail some special cases of this dependence for 
a— \ (solid lines), a — 1.25 (dashed lines) and a — 0.75 
(dash-dotted fines). 

We note that, in terms of the general equation J^l, the 
cases of a < 1 and a > 1 are tantamount to each other, 
as one may divide the equation by a, mutually rename 
the vertical and horizontal indices (n and m), and then 
rescale the equation to the form with a replaced by l/a. 
However, this transformation is not possible once we fix 
A = 1, which is why we report results below for both 
a > 1 and a < 1. 

For a = 1, an eigenvalue bifurcates from the edge of 
the continuous spectrum at e w 0.445, and with further 
increase of e it moves towards the origin of the spec- 
tral plane (Ar,Ai) (the subscripts denote the real and 
imaginary part of the eigenvalue) . It becomes unstable, 
reaching the origin at e « 1.006. For a = 1.25, the first 
bifurcation occurs at e ~ 0.398, and the instability sets in 
at e « 0.896, whereas for a — 0.75 the respective critical 
points (the appearance of the eigenvalue, and its passage 
into the instability region) are found at e « 0.511 and 
1.156, respectively. Notice that these results are quite 
natural since, as a — > 0, the system becomes nearly one- 
dimensional, hence we expect the destabilization point to 
approach its ID counterpart. Thus, as the ID discrete 
solitons are well-known to be stable up to the continuum 
limit, one may expect that — > oo for a 0. The bot- 
tom panel of Fig.|2]shows an example of a discrete soliton 
for a = 1.5 and e = 1. Although the anisotropy is hardly 
observed in this case, it can be traced nevertheless; in 
particular, uifl — 0.785, and mo,i — 0.579. 

Similar results can be obtained for on-site vortices 
(discrete vortex solitons) with the topological charge 
S = 1. In this section, we consider the solitons in 
the form of the so-called "vortex cross", with ui^ = 1, 
1*0,1 = exp(i7r/2) = i, u_i,o — exp(i7r) = —1, uo,-i = 
exp(z37r/2) = —i (and Unn = 0, at the central point), 
excited in the AC limit There are interesting vari- 
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FIG. 2: The soliton line in the top panel shows the critical 
value of e (the border between stable and unstable discrete 
solitons) as a function of a; the dashed line is e = The 
middle panel shows how the real and imaginary parts of the 
stability eigenvalue, \r and A^, depend on e for ct — 1.25, 1, 
and 0.75 (dashed, solid, and dash-dotted curves, respectively). 
The bottom panel shows an example of the discrete soliton 
found for e — 1 and a = 1.5. 



ations to this problem, in comparison with the funda- 
mental soliton. In particular, the respective instability 
mechanism is different, as it is caused by an eigenvalue 
bifurcating from the origin in the spectral plane for e 7^ 0, 
and eventually (upon parametric continuation) colliding 
with the edge of the continuous spectrum (or an eigen- 
value bifurcating from the continuous spectrum). The 
collision gives rise to a quartet of eigenvalues, through 
the so-called Hamiltonian-Hopf bifurcation [2J|. In the 
isotropic case (a = 1), it is known that this instability 
sets in at ~ 0.39 while, in the present anisotropic 
model, we have found that £„■ ~ 0.325 for a = 1.3 and 



Ecr ~ 0.429 for a = 0.7. The respective two-parameter 
diagram (£„•, a) is shown in the top panel of Fig. 13 The 
cases of a = 1.3, 1, and 0.7 (dashed, solid and dash- 
dotted, dashed-dotted curves, respectively) are shown in 
the middle panel. The bottom panel of the figure il- 
lustrates the squared-amplitude profile of the discrete 
vortex for a = 0.2 and e — 0.5. The sites (1,0) and 
(0,1) have the squared amplitudes |ui_oP = 1.934 and 
|uo,ip = 2.057, respectively. Notice also that as a ^ 0, 
a quasi- ID situation is again approached, where the so- 
called twisted-localized mode (TLM) 25j configuration 
(alias an odd soliton) is a counterpart of the 2D vortex. 
As one would expect, the critical point of the instability 
departs from the value e'^^^ ~ 0.39, corresponding to 
the isotropic 2D case, towards the value corresponding 
to the stability border of the ID TLM solitons, which is 



0.433. 



IV. FUNDAMENTAL VORTEX SQUARES 

For the discrete solitons examined so far, the differ- 
ence between the isotropic and non-isotropic cases has 
not been particularly dramatic; the anisotropy chiefly 
entailed a smooth deformation of the instability-onset 
scenarios known for the isotropic case. Therefore, the 
dynamical evolution triggered by the instability is natu- 
rally expected to be similar to that in previously studied 
isotropic cases iZi ^ ^, . 

Now we will give an example where the instability sce- 
nario and dynamics are very different from their isotropic 
counterparts. We focus, in particular, on the off site- 
centered vortex (alias "vortex square" ) The vortex- 
square contours are characterized by their size M , which 
is the number of lattice bonds that each side of the square 
contour contains in the AC-limit pattern, from which the 
solution family stems. Hence, the vortex square based, in 
the AC limit, on the set of sites (0, 0), (1, 0), (1, 1), (0, 1) 
is the M = 1 contour. The configuration with 5 = 1 is 
written on this set by lending the four sites the phases 
0, 7r/2, TT and 37r/2, respectively. The persistence of such 
configurations, as was discussed in detail in Ref. [2^ . 
is determined by whether secular conditions (obtained 
from the Lyapunov-Schmidt theory p^*). excluding the 
projection of eigenvectors in the kernel of the lineariza- 
tion at e = to the solution at finite e, are satisfied. In 
the isotropic case, to leading order (0(e)), these secular 
conditions are found to be 



= f{ei) 



9i - di+i) + sinidi 



n-1 



(21) 



for I = 1,...,N (with periodic boundary conditions), 
where N = 4Af is the number of sites participating in 
the contour and 9i are their respective phases [cf. Eqs. 
(3.1)~(3.2) of Ref. jH]. 

One can then apply similar arguments to the present 
setting and derive modified persistence criteria for the 
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FIG. 3: The top panel shows the critical value of e separat- 
ing the stable and unstable discrete vortices (on-site-centered 
ones, alias vortex crosses) with 5* = 1 as a function of a. 
The middle panel shows how the real and imaginary parts 
of the eigenvalue leading to the instability depend on e for 
a = 1.3, 1, and 0.7 (dashed, solid and dash-dotted curves, 
respectively). Notice that, for a = 1.3, there is a secondary 
instability arising for e > 0.433. The bottom panel shows 
the squared-absolute-value profile of the discrete vortex for 
e = 0.5 and a = 0.2. 



anisotropic model. For M = 1, they are 



= f{Oi) = 



asm[tli - Ui+i) + sin(«; - yj-ij 

I = 2fc + l,fc = 0, 1 

sin(6'; - 6';+i) + asm{9i ~ 9i^i) 

I = 2k, k = 1,2. 



(22) 

While Eqs. H22|l may seem a moderate modification of 
()21|l . there is a crucial (for stability purposes) difference. 
Indeed, consider the linearization around the S = 1 so- 



lution according to Eq. It was proved in Ref. |2g 
that the Jacobian matrix of the reduced set of Eqs. H22|l. 
defined through Jik = dfi/d9k, determines leading-order 
corrections to A^ — 1 eigenvalue pairs bifurcating from the 
origin [one pair stays at the origin due to the invariance 
of Eq. Q with respect to the phase shift], since these 
eigenvalues satisfy the equation: 



(23) 



with /i; the corresponding eigenvalues of the reduced 
N X N Jacobian Jik- It is further easy to check that, 
for the vortex square with S = 1 and M = 1, the en- 
tire Jacobian matrix consists of zeros. More generally, 
as shown in Ref. |2^, this is the case for the square 
vortices of size M with charge S — M, which for that 
reason were termed "super-symmetric" vortices. Obvi- 
ously, to determine the stability of the vortices in this 
special case, one needs to go to higher-order expansions. 
Typically, second-order reductions will yield a non-trivial 
result for the stability of such super-symmetric configu- 
rations, leading to eigenvalue dependences A; oc e [rather 
than Xj oc ^/e, as dictated by Eq. H23|l in the generic 
case] . 

The key variation to this theme stemming from the 
presence of the anisotropy is that the matrix Jj^ has 
generically non- vanishing elements in the lowest approxi- 
mation for a 7^ 1; in other words, the isotropic lattice is a 
degenerate one for the supersymmetric solitons, and arbi- 
trarily weak anisotropy lifts this degeneracy. As a result, 
the eigenvalue bifurcations occur, typically, at the lead- 
ing order, rather than at the second-order perturbation 
expansions, which was the case in the isotropic model. 
More strikingly, considering a specific example, such as 
for a = 1.05 (a very weak deviation from the isotropic 
case), we find that the relevant angles (in radians) satis- 
fying the conditions ^ are 6*1 = -0.0229, 6*2 = 1.8577, 
6*3 = 3.4285, and 64 = 4.6895; the corresponding 4 x 4 Ja- 
cobian has two zero eigenvalues (one of which will split to 
order 0(e), see below) and two nonzero ones, ±0.6403. 
From the existence of the positive eigenvalue and from 
Eq. it immediately follows that the S = M = 1 

configuration is immediately unstable (for all values of 
e) . This is in complete contrast with the super-symmetric 
vortex in the isotropic model, which has two imaginary 
eigenvalue pairs (bifurcating at the second-order reduc- 
tion), A « ±2ie, and is linearly stable for e < Ec ~ 0.38. 

From here, we conclude that the anisotropy can play a 
critical role in destabilizing configurations that would be 
very robust ones in the isotropic limit. Furthermore, this 
can happen arbitrarily close to the isotropic limit (that 
turns out to be a very delicate one), given the nature 
of the argument presented above. We also note in pass- 
ing that in the anisotropic case examined above, there 
is yet another real eigenvalue pair which is A ~ ±3e for 
small e (this pair stems from the higher-order reduction, 
in agreement with the prediction of the reduced Jaco- 
bian). These two eigenvalue pairs eventually collide at 
e = 0.057, resulting in a Hamiltonian Hopf bifurcation 
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0.08 0.1 



FIG. 4; For the 5 = 1 super-symmetric square vortex (one 
with the size M = 1) , two real stability eigenvalues are shown 
as functions of e for a = 1.05. The numerical and analytical 
results (see text) are displayed, respectively, by the solid lines 
and dashed lines. 



to an eigenvalue quartet which is present in the stability 
spectrum at e > 0.057. This phenomenology is shown in 
Fig. 21 The leading-order prediction for the most unsta- 
ble eigenvalue is in good agreement with the full numeri- 
cal result for small values of e. For higher values of e, the 
second-order corrections that we do not examine here in 
detail come into play and lead to the Hamiltonian Hopf 
bifurcation. 

To directly compare the dynamics between the 
isotropic and weakly anisotropic (yet unstable) case for 
the super-symmetric vortex, we have performed numeri- 
cal simulations. Detailed simulations are reported in this 
work only for the super-symmetric cases (see also the 
next section), since for all other states anisotropy oper- 
ates as a regular perturbation, see above; as a result, 
instabilities of the other states may be shifted due to the 
anisotropy, but structurally the phenomenology remains 
the same. 

For the delicate super-symmetric vortex square, the 
dynamics altered by the anisotropy is indeed found to 
be dramatically different from the isotropic case. This 
is illustrated by Fig. [Sj for the vortex square with S = 
M = 1, carried (in the AC limit) by 4 sites. The time 
dynamics of the squared absolute value of the field at the 
main sites is shown in the figure for a weakly anisotropic 
model, with a = 1.05, and its isotropic counterpart (top 
and bottom panels, respectively). Stark contrast between 
the instability developing for i > 50 in the former case, 
versus the complete stability for all times in the latter 
(isotropic) system, is obvious (notice the difference in 
the scales of vertical axes between the two panels). In 
the linear approximation, these results are well predicted 
by the above theory. 
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FIG. 5: The dynamics of an initially very weakly perturbed 
super-symmetric vortex with S = M = 1, principally based 
on four lattice sites that form an elementary cell (the sites 
are labeled as 1,2,3,4). The time evolution of the squared 
absolute value of the fields at these sites is shown in the top 
panel for a weakly anisotropic model, with a — 1.05, and for 
its isotropic counterpart (a — 1) in the bottom panel. In 
both cases, the same uniformly distributed, random initial 
perturbation of amplitude 10~* was added to the solution at 
t = to excite possible instabilities. Clearly, the vortex on 
the weakly anisotropic lattice becomes unstable at t > 50, 
while in the isotropic case the perturbation remains bounded 
and small at all times. In these examples, the intersite lattice 
coupling constant is e = 0.025. 



V. HIGHER-ORDER VORTICES 

We now give a summary of results for vortices with 
higher values of the topological charge. First, we con- 
sider the S ^ M = 2 super-symmetric vortex populating 
the sites (1,0), (1,1), (0,1), (-1,1), (-1,0), (-1,-1), 
(0, —1), (1, —1) in the AC limit, with a phase shift of 7r/2 
between adjacent sites (in the isotropic model). The lat- 
ter provides for a total phase gain of An around a closed 
path surrounding the origin. This type of the configura- 
tion with S = M = 2 was identified in Ref. "2^ as pos- 
sessing a real eigenvalue pair with Xr — ± m 
excellent agreement with numerical computations. How- 
ever, the presence of the small anisotropy for a 7^ 1 again 
strongly affects the vortex for reasons similar to the ones 
presented above. In this case, the reductions leading to 
the perturbed dynamics in the anisotropic model are de- 
scribed by the following persistence conditions: 

sin(6l, - 6l,+i) + sin(6'; - 6i/_i) 

i = 2fc + l,/c = 0,1,2,3, 



= f{Ol) EE { 



asin{9i - 0i+i) + sin(6'/ - 0;_i) 

i = 4fc + 2,A: = 0, 1, 

sin(6'; - 6*;+!) + asin(6'; - 0(_i) 

i = 4fc-|-4,/c = 0, 1, 



(24) 
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0.18 




FIG. 6: For the S — M = 2 supersymmetric vortex, the three 
real eigenvalues are displayed as functions of e for a = 1.05. 
The solid and dashed lines depict the numerical and analytical 
results. 



cf. Eqs. (|22|l . In this expression 9i is the phase of the 
field at each of the eight above-mentioned sites (where, 
in the order the sites were mentioned, the correspond- 
ing index is / — 1,2, ... ,8). Furthermore, as discussed 
above, the analysis performed in Ref. [2^ can be used to 
show that the hnear stability eigenvalues for such a vor- 
tex soliton will be given, to the leading order, by Eq. (|23|l . 
Using this prediction, even in the weakly anisotropic 
case (e.g., for a = 1.05) one finds that the correspond- 
ing 8x8 Jacobian possesses three real 0{y/£) eigenval- 
ues, which result in an instability (contrary to the single 
real 0{e) eigenvalue in the a = 1 case). Hence, once 
again, the anisotropy results in a significant destabiliza- 
tion of the super-symmetric vortex, in comparison to the 
isotropic model. As a specific example, we show in Fig. 
El the situation for a = 1.05. The solution of Eqs. if^ 
yields Oi = 0.218, 6*2 = 1.967, 6*3 = 3.182, 64 = 4.397, 
6*5 = 6.145, 9e = 7.894, 6*7 9.109, 6*8 = 11.036, which, 
in turn, results in a Jacobian with the 3 real eigenval- 
ues n = {1.0145,0.5357,0.2391}. The comparison of the 
numerical prediction for the eigenvalue dependence on e 
versus the corresponding analytical prediction (solid and 
dashed lines, respectively) based on the above results is 
given in Fig. demonstrating a very good agreement 
between the two. 

To highlight the substantial differences between the dy- 
namics in the isotropic and anisotropic models, we have 
performed numerical simulations of the super-symmetric 
vortex with S = M ~ 2. In this case, the evolution of the 
field at the eight basic sites is shown in the top panel of 
Fig.[7|for a = 1.05, and in the bottom panel for a = 1. In 
the former case, for the coupling strength e = 0.015 con- 
sidered here, the three unstable eigenvalues for a — 1.05 
are A = 0.1688, A 0.1258 and A = 0.0855, while in the 
latter case (isotropic model) , the only unstable eigenvalue 
is a much smaller one, A = 0.0146. Naturally, we observe 
the instability setting in much earlier in the anisotropic 
model (at t ^ 30) than in the isotropic one (at t ^ 160). 

One may be wondering whether the strong dynamical 




100 200 , 300 400 



FIG. 7: Same as in Fig. |S| but for the supersymmetric vor- 
tex of the S — M = 2 type. The four lines (solid, dashed, 
dotted, and dash-dotted) and four symbols (circles, pluses, 
stars and triangles) are used to denote the squared absolute 
values of the field at the eight sites carrying the vortex in the 
anisotropic (top) model and its isotropic counterpart (bot- 
tom) for e = 0.015. 

effect of the weak anisotropy should be attributed to the 
super-symmetry of the vortex, or maybe just the specific 
type of contour which carries the vortex. To check this, 
we have also considered the vortex with 5* = 3 sitting on 
the same M = 2 contour. Given the lack of the super- 
symmetry in the latter case, the bifurcation of the rele- 
vant 7 (= — 1) eigenvalue pairs occurs at the leading- 
order reduction and all of them are proportional to 
More specifically, for the largest pair in the isotropic case 
(for instance), the proportionality factor is 2.3784. In the 
anisotropic case with a = 1.05, the seven pairs remain 
on the imaginary axis, being slightly perturbed due to 
a 1. For instance, the largest one among them is now 
A = ±2.3943 • i,/e. On the other hand, for a = 0.95, the 
largest eigenvalue pair is A = ±2.3647 • iy^. This also is 
in line with our above results on the fundamental discrete 
soliton and vortex cross, since it indicates that, for a > 1, 
the collision of this eigenvalue with the continuous spec- 
trum (which leads to the Hamiltonian Hopf bifurcation) 
will occur at smaller e, the opposite being true for a < 1. 
Hence the stability diagram of the S ^ 3, M = 2 vortex 
square is quite similar to that shown for the fundamental 
soliton and vortex cross in Figs. El and |2I (therefore, it is 
not shown here). 



VI. CONCLUSIONS 

In this work, we have examined effects of anisotropy 
on lattice nonlinear dynamical systems supporting dis- 
crete solitons and vortices. The two-dimensional dis- 
crete nonlinear Schrodinger equation was used as a 
paradigm model. The variational approximation was 
developed for fundamental solitons, showing (by means 
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of the Vakhitov-Kolokolov criterion) that broad quasi- 
continuum ones are unstable, while strongly anisotropic 
solitons are stable. By means of numerical methods, we 
have found that usual localized states, such as the fun- 
damental discrete solitons and vortex crosses, are only 
mildly affected by the anisotropy, which results in a modi- 
fied stability region (reduced when one direction features 
a stronger coupling than the isotropic limit, and aug- 
mented when the coupling along this direction is weaker). 
General phenomenology for such states is similar to that 
for their counterparts on the isotropic lattice. 

The main finding reported in the present work is that 
the assumption about mild deformation of the stability 
region induced by weak anisotropy is not valid for the del- 
icate super-symmetric vortex states residing on square 
contours, in the case when the vorticity S is equal to 
the contour's size M. In this special case, the degener- 
acy of the leading-order existence conditions (dictated by 
Lyapunov-Schmidt theory) specific to the isotropic case is 
broken by the anisotropy. This, in turn, results in a dra- 
matically different behavior (as a function of the intersite 
coupling constant) of the corresponding linear stability 
eigenvalues, in terms of both the order of their bifur- 
cation, and the number of real eigenvalues. As a conse- 
quence, the supersymmetric vortex-square structure that 
was marginally stable in the isotropic case is found to be 
strongly unstable even on the weakly anisotropic lattice. 
Similarly, the supersymmetric vortex with 5 = A/ = 2 is 
found to be much more unstable in the anisotropic case 



in comparison to its isotropic counterpart. 

The most natural systems for experimental observa- 
tion of the results predicted in this work are deep optical 
lattices trapping BECs, and bundled sets of nonlinear 
optical waveguides (the latter have been recently created 
experimentally [l^). Anisotropic lattices can also be in- 
duced in photorefractive media, but this medium should 
be considered separately, in view of the different (sat- 
urable) character of the optical nonlinearity in this case. 
Such investigations are currently in progress and will be 
reported elsewhere. 

A further natural extension of this work would be 
to examine effects of anisotropy in three-dimensional 
lattices on discrete solitons, vortices, dipoles and 
quadrupoles of various types, octupoles, and more exotic 
localized configurations, that were recently investigated 
for the isotropic case in Refs. p^. 
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